*Import the sample data: 
use "$root/Data/Produced/EV_sample.dta", clear

*Generate the choice data (i.e. each household had each available car at disposal in a given year)

keep PID car_combo year_reg

fillin PID car_combo

gen choice=(_fillin==0)

bys PID (_fillin): replace year_reg=year_reg[_n-1] if missing(year_reg)

drop _fillin

merge m:1 car_combo year_reg using "$root/Data/Produced/sample_car_chars.dta", keepusing(MAKName drivetype car_size car_height et_verbrauch el_verbrauch env_cat PRHnewprice weight TYPKW car_type_FE* )

drop _merge

merge m:1 PID using "$root/Data/Produced/sample_soc_chars.dta", keepusing(average_use_var TOT_EINK TOTVERM PVhh homeowner hhsize age_group range5 dist_to_EV urban3)

drop _merge

*Now we also add the instrumental variables: 
merge m:1 car_combo year_reg using "$root/Data/Produced/diff_instruments.dta"

drop _merge

*Now we have the data ready to conduct choice estimation.
 
*We first produce all the required additional variables: 

*Generate the dummies for the drivetype: 

gen diesel=(drivetype==2)
gen ev=(drivetype==3)
gen hyrid=(drivetype==4)

*Generate the price variable which is measured in 1000: 
gen price=PRHnewprice/1000

*rename the height variable: 
rename car_height height

*rename the KW variable: 
rename TYPKW typkw

*Generate the variable average_use which is assumed as 12000km for all users and 16000km for diesel users: 
gen average_use=12000 
replace average_use=16000 if diesel==1

*Generate the brand_FE for the top10 brands: 
gen brand_10=MAKName
replace brand_10="other US" if MAKName=="cadillac" |MAKName=="chevrolet"  | MAKName=="jeep" |MAKName=="tesla"
replace brand_10="other asia" if (MAKName=="honda") |(MAKName=="hyundai") |(MAKName=="infiniti") |(MAKName=="kia") | (MAKName=="lexus") |(MAKName=="mazda") |(MAKName=="mitsubishi") |(MAKName=="nissan") |(MAKName=="ssangyong") |(MAKName=="subaru")  |(MAKName=="toyota")
 
replace brand_10="other german" if (MAKName=="audi") |(MAKName=="porsche") |(MAKName=="smart")   

replace brand_10="other french" if (MAKName=="citroen") | (MAKName=="ds") 

replace brand_10="other Italian" if (MAKName=="alfa-romeo") |(MAKName=="fiat") |(MAKName=="maserati") 
replace brand_10="other EU" if (MAKName=="dacia") |(MAKName=="jaguar") |(MAKName=="land-rover") |(MAKName=="mini") |(MAKName=="volvo") 

tab brand_10, gen(brand_FE)

drop brand_10 MAKName brand_FE1


*Generate the wealth groups and interaction terms with wealth: 
xtile wealth_group=TOTVERM, nquantiles(4)

tab wealth_group, g(wealth_quart)

gen ev_wealth2=ev*wealth_quart2
gen ev_wealth3=ev*wealth_quart3
gen ev_wealth4=ev*wealth_quart4

gen price_wealth2=price*wealth_quart2
gen pice_wealth3=price*wealth_quart3
gen price_wealth4=price*wealth_quart4

gen wealth=TOTVERM/1000

drop wealth_group1 TOTVERM

*Generate the interaction terms between EV and urbanity: 
gen ev_urban2=(ev==1)&(urban3==2)
gen ev_urban3)(ev==1)&(urban3==3)

*Generate the interaction terms between EV and other socio-economic characteristics: 
gen ev_dist_ev=ev*dist_to_EV
gen ev_range5=ev*range5
gen ev_year2=(ev==1)&(year_reg==2018)
gen ev_year3=(ev==1)&(year_reg==2019) 
gen ev_pv=ev*PVhh
gen ev_ho=ev*homeowner

*Drop the no longer needed vars: 
drop dist_to_EV range5 PVhh homeowner

**Generate the size & household size interactions: 
tab hhsize, g(hh)
gen size_hh2=car_size*hh2 
gen size_hh3=car_size*hh3 
gen size_hh4=car_size*hh4
gen size_hh5=car_size*hh5

drop hh1 hh2 hh3 hh4 hh5 hhsize
*Generate the power age_group interactions: 
tab age_group, g(age)
gen age_kw2=age2*typkw
gen age_kw3=age3*typkw

drop age1 age2 age3 age_group

*Generate the variable environmentally friendly: 
gen env_friendly=(env_cat>=3)

** Calculate the variable drive cost, which measures costs per 100km: 
**Generate drivecosts for 100km
gen drive_cost=.

*There was a comma missing in el_verbrauch
replace el_verbrauch=el_verbrauch/10

*For gasoline cars: 1.51 (2017), 1.63 (2018), 1.60 (2019)
replace drive_cost=et_verbrauch*1.51 if drivetype==1 & year_reg<=2017
replace drive_cost=et_verbrauch*1.63 if drivetype==1 & year_reg==2018
replace drive_cost=et_verbrauch*1.60 if drivetype==1 & year_reg==2019

*For diesel cars: 1.58 (2017), 1.74 (2018), 1.74 (2019)
replace drive_cost=et_verbrauch*1.58 if drivetype==2 & year_reg<=2017
replace drive_cost=et_verbrauch*1.74 if drivetype==2 & year_reg==2018
replace drive_cost=et_verbrauch*1.74 if drivetype==2 & year_reg==2019

*For electronic cars 0.2225 (2017), 0.2272 (2018), 0.2277 (2019):
replace drive_cost=el_verbrauch*0.2225 if drivetype==3 & year_reg<=2017
replace drive_cost=el_verbrauch*0.2272 if drivetype==3 & year_reg==2018
replace drive_cost=el_verbrauch*0.2277 if drivetype==3 & year_reg==2019

*For hybrid car gasoline and electric assumed
replace drive_cost=et_verbrauch*1.51 + el_verbrauch*0.2225 if drivetype==4 & year_reg<=2017
replace drive_cost=et_verbrauch*1.63 + el_verbrauch*0.2272 if drivetype==4 & year_reg==2018
replace drive_cost=et_verbrauch*1.60 + el_verbrauch*0.2277 if drivetype==4 & year_reg==2019

**Add the program to compute the vehicle tax for one year: 

*Program the tax for one year: 
gen car_tax=.
replace car_tax=240*(weight/1000) if weight<1000
replace car_tax=240+0.2064*(weight-1000) if 1000<=weight<=2000&drivetype!=3
replace car_tax=240+206.4+0.1775*(weight-2000) if 2000<weight<=3000&drivetype!=3
replace car_tax=240+206.4+177.5+0.1527*(weight-3000) if weight>3000&drivetype!=3

*Electric vehicles pay only half the tax rate (120.- per kg)
replace car_tax=120*(weight/1000) if weight<1000
replace car_tax=120+0.1032*(weight-1000) if 1000<=weight<=2000&drivetype==3
replace car_tax=120+103.2+0.08875*(weight-2000) if 2000<weight<=3000&drivetype==3
replace car_tax=120+103.2+88.8+0.07633*(weight-3000) if weight>3000&drivetype==3

*In the first 3 years of ownership, Cat A & B cars as well as electric cars are taxed at a reduced rate
replace car_tax=car_tax*0.6 if drivetype!=3 & env_cat==1 
replace car_tax=car_tax*0.8 if drivetype!=3 & env_cat==2
replace car_tax=car_tax*0.4 if drivetype==3 

*Calculate the annual drive cost for both variable and non-variable kilometres: 
gen drive_cost_ann=average_use*drive_cost/100
gen drive_cost_ann_variable=average_use_var*drive_cost/100


**Now we calculate the different net present values of the variable costs and taxes using different assumptions: 
local r=0.06


gen pv_drive_cost=drive_cost_ann+drive_cost_ann*((1-((1+`r')^(-(14))))/`r') 




**Assuming a car-lifetime of 15 years and an interest rate of 6% as the literature does we calculate the pv of the tax revenue that is still owed (as of 2019): 
local r=0.06

gen pv_tax=.

replace pv_tax=car_tax+(car_tax)/((1+`r')^1)+(car_tax)/((1+`r')^2)+(car_tax)/((1+`r')^3)+(car_tax)/((1+`r')^(4)) +(car_tax)/((1+`r')^(5)) +(car_tax)/((1+`r')^(6)) +(car_tax)/((1+`r')^(7)) +(car_tax)/((1+`r')^(8)) +(car_tax)/((1+`r')^(9)) +(car_tax)/((1+`r')^(10)) +(car_tax)/((1+`r')^(11)) +(car_tax)/((1+`r')^(12))+(car_tax)/((1+`r')^(13)) +(car_tax)/((1+`r')^(14)) if drivetype!=3 & env_cat>2

replace pv_tax=car_tax+(car_tax)/((1+`r')^1)+(car_tax)/((1+`r')^2)+(car_tax)/((1+`r')^3)+(car_tax/0.4)/((1+`r')^(4)) +(car_tax/0.4)/((1+`r')^(5)) +(car_tax/0.4)/((1+`r')^(6)) +(car_tax/0.4)/((1+`r')^(7)) +(car_tax/0.4)/((1+`r')^(8)) +(car_tax/0.4)/((1+`r')^(9)) +(car_tax/0.4)/((1+`r')^(10)) +(car_tax/0.4)/((1+`r')^(11)) +(car_tax/0.4)/((1+`r')^(12))+(car_tax/0.4)/((1+`r')^(13)) +(car_tax/0.4)/((1+`r')^(14)) if drivetype==3

replace pv_tax=car_tax+(car_tax)/((1+`r')^1)+(car_tax)/((1+`r')^2)+(car_tax)/((1+`r')^3)+(car_tax/0.6)/((1+`r')^(4)) +(car_tax/0.6)/((1+`r')^(5)) +(car_tax/0.6)/((1+`r')^(6)) +(car_tax/0.6)/((1+`r')^(7)) +(car_tax/0.6)/((1+`r')^(8)) +(car_tax/0.6)/((1+`r')^(9)) +(car_tax/0.6)/((1+`r')^(10)) +(car_tax/0.6)/((1+`r')^(11)) +(car_tax/0.6)/((1+`r')^(12))+(car_tax/0.6)/((1+`r')^(13)) +(car_tax/0.6)/((1+`r')^(14)) if drivetype!=3 & env_cat==1

replace pv_tax=car_tax+(car_tax)/((1+`r')^1)+(car_tax)/((1+`r')^2)+(car_tax)/((1+`r')^3)+(car_tax/0.8)/((1+`r')^(4)) +(car_tax/0.8)/((1+`r')^(5)) +(car_tax/0.8)/((1+`r')^(6)) +(car_tax/0.8)/((1+`r')^(7)) +(car_tax/0.8)/((1+`r')^(8)) +(car_tax/0.8)/((1+`r')^(9)) +(car_tax/0.8)/((1+`r')^(10)) +(car_tax/0.8)/((1+`r')^(11)) +(car_tax/0.8)/((1+`r')^(12))+(car_tax/0.8)/((1+`r')^(13)) +(car_tax/0.8)/((1+`r')^(14)) if drivetype!=3 & env_cat==2

*Also generate the pv of driving costs: 

gen var_cost=(pv_tax+pv_drive_cost)/1000

*Generate the drive costs with variable amount of KMs: 
*Generate the driving costs with constant amount of kilometres: 
**Generate the constant costs: 
local r=0.06


gen pv_drive_cost_var=drive_cost_ann_variable+drive_cost_ann_variable*((1-((1+`r')^(-(14))))/`r') 


gen var_cost_variable=(pv_tax+pv_drive_cost_var)/1000


**Assuming a car-lifetime of 15 years and an interest rate of 2% for the sensitivity check low:
local r=0.02

gen pv_tax_low=.

replace pv_tax_low=car_tax+(car_tax)/((1+`r')^1)+(car_tax)/((1+`r')^2)+(car_tax)/((1+`r')^3)+(car_tax)/((1+`r')^(4)) +(car_tax)/((1+`r')^(5)) +(car_tax)/((1+`r')^(6)) +(car_tax)/((1+`r')^(7)) +(car_tax)/((1+`r')^(8)) +(car_tax)/((1+`r')^(9)) +(car_tax)/((1+`r')^(10)) +(car_tax)/((1+`r')^(11)) +(car_tax)/((1+`r')^(12))+(car_tax)/((1+`r')^(13)) +(car_tax)/((1+`r')^(14)) if drivetype!=3 & env_cat>2

replace pv_tax_low=car_tax+(car_tax)/((1+`r')^1)+(car_tax)/((1+`r')^2)+(car_tax)/((1+`r')^3)+(car_tax/0.4)/((1+`r')^(4)) +(car_tax/0.4)/((1+`r')^(5)) +(car_tax/0.4)/((1+`r')^(6)) +(car_tax/0.4)/((1+`r')^(7)) +(car_tax/0.4)/((1+`r')^(8)) +(car_tax/0.4)/((1+`r')^(9)) +(car_tax/0.4)/((1+`r')^(10)) +(car_tax/0.4)/((1+`r')^(11)) +(car_tax/0.4)/((1+`r')^(12))+(car_tax/0.4)/((1+`r')^(13)) +(car_tax/0.4)/((1+`r')^(14)) if drivetype==3

replace pv_tax_low=car_tax+(car_tax)/((1+`r')^1)+(car_tax)/((1+`r')^2)+(car_tax)/((1+`r')^3)+(car_tax/0.6)/((1+`r')^(4)) +(car_tax/0.6)/((1+`r')^(5)) +(car_tax/0.6)/((1+`r')^(6)) +(car_tax/0.6)/((1+`r')^(7)) +(car_tax/0.6)/((1+`r')^(8)) +(car_tax/0.6)/((1+`r')^(9)) +(car_tax/0.6)/((1+`r')^(10)) +(car_tax/0.6)/((1+`r')^(11)) +(car_tax/0.6)/((1+`r')^(12))+(car_tax/0.6)/((1+`r')^(13)) +(car_tax/0.6)/((1+`r')^(14)) if drivetype!=3 & env_cat==1

replace pv_tax_low=car_tax+(car_tax)/((1+`r')^1)+(car_tax)/((1+`r')^2)+(car_tax)/((1+`r')^3)+(car_tax/0.8)/((1+`r')^(4)) +(car_tax/0.8)/((1+`r')^(5)) +(car_tax/0.8)/((1+`r')^(6)) +(car_tax/0.8)/((1+`r')^(7)) +(car_tax/0.8)/((1+`r')^(8)) +(car_tax/0.8)/((1+`r')^(9)) +(car_tax/0.8)/((1+`r')^(10)) +(car_tax/0.8)/((1+`r')^(11)) +(car_tax/0.8)/((1+`r')^(12))+(car_tax/0.8)/((1+`r')^(13)) +(car_tax/0.8)/((1+`r')^(14)) if drivetype!=3 & env_cat==2

*Also generate the pv of driving costs: 
gen pv_drive_cost_low=drive_cost_ann + drive_cost_ann*((1-((1+`r')^(-(14))))/`r') 

gen var_cost_low=(pv_tax_low+pv_drive_cost_low)/1000

drop pv_tax_low pv_drive_cost_low

**Assuming a car-lifetime of 15 years and an interest rate of 10% for the sensitivity check high:
local r=0.1

gen pv_tax_high=.

replace pv_tax_high=car_tax+(car_tax)/((1+`r')^1)+(car_tax)/((1+`r')^2)+(car_tax)/((1+`r')^3)+(car_tax)/((1+`r')^(4)) +(car_tax)/((1+`r')^(5)) +(car_tax)/((1+`r')^(6)) +(car_tax)/((1+`r')^(7)) +(car_tax)/((1+`r')^(8)) +(car_tax)/((1+`r')^(9)) +(car_tax)/((1+`r')^(10)) +(car_tax)/((1+`r')^(11)) +(car_tax)/((1+`r')^(12))+(car_tax)/((1+`r')^(13)) +(car_tax)/((1+`r')^(14)) if drivetype!=3 & env_cat>2

replace pv_tax_high=car_tax+(car_tax)/((1+`r')^1)+(car_tax)/((1+`r')^2)+(car_tax)/((1+`r')^3)+(car_tax/0.4)/((1+`r')^(4)) +(car_tax/0.4)/((1+`r')^(5)) +(car_tax/0.4)/((1+`r')^(6)) +(car_tax/0.4)/((1+`r')^(7)) +(car_tax/0.4)/((1+`r')^(8)) +(car_tax/0.4)/((1+`r')^(9)) +(car_tax/0.4)/((1+`r')^(10)) +(car_tax/0.4)/((1+`r')^(11)) +(car_tax/0.4)/((1+`r')^(12))+(car_tax/0.4)/((1+`r')^(13)) +(car_tax/0.4)/((1+`r')^(14)) if drivetype==3

replace pv_tax_high=car_tax+(car_tax)/((1+`r')^1)+(car_tax)/((1+`r')^2)+(car_tax)/((1+`r')^3)+(car_tax/0.6)/((1+`r')^(4)) +(car_tax/0.6)/((1+`r')^(5)) +(car_tax/0.6)/((1+`r')^(6)) +(car_tax/0.6)/((1+`r')^(7)) +(car_tax/0.6)/((1+`r')^(8)) +(car_tax/0.6)/((1+`r')^(9)) +(car_tax/0.6)/((1+`r')^(10)) +(car_tax/0.6)/((1+`r')^(11)) +(car_tax/0.6)/((1+`r')^(12))+(car_tax/0.6)/((1+`r')^(13)) +(car_tax/0.6)/((1+`r')^(14)) if drivetype!=3 & env_cat==1

replace pv_tax_high=car_tax+(car_tax)/((1+`r')^1)+(car_tax)/((1+`r')^2)+(car_tax)/((1+`r')^3)+(car_tax/0.8)/((1+`r')^(4)) +(car_tax/0.8)/((1+`r')^(5)) +(car_tax/0.8)/((1+`r')^(6)) +(car_tax/0.8)/((1+`r')^(7)) +(car_tax/0.8)/((1+`r')^(8)) +(car_tax/0.8)/((1+`r')^(9)) +(car_tax/0.8)/((1+`r')^(10)) +(car_tax/0.8)/((1+`r')^(11)) +(car_tax/0.8)/((1+`r')^(12))+(car_tax/0.8)/((1+`r')^(13)) +(car_tax/0.8)/((1+`r')^(14)) if drivetype!=3 & env_cat==2
*Also generate the pv of driving costs: 
gen pv_drive_cost_high=drive_cost_ann + drive_cost_ann*((1-((1+`r')^(-(14))))/`r') 

gen var_cost_high=(pv_tax_high+pv_drive_cost_high)/1000

drop pv_tax_high pv_drive_cost_high


**Assuming a car-lifetime of 25 years and an interest rate of 6% for the sensitivity check long:
local r=0.06

gen pv_tax_long=.

replace pv_tax_long=car_tax+(car_tax)/((1+`r')^1)+(car_tax)/((1+`r')^2)+(car_tax)/((1+`r')^3)+(car_tax)/((1+`r')^(4)) +(car_tax)/((1+`r')^(5)) +(car_tax)/((1+`r')^(6)) +(car_tax)/((1+`r')^(7)) +(car_tax)/((1+`r')^(8)) +(car_tax)/((1+`r')^(9)) +(car_tax)/((1+`r')^(10)) +(car_tax)/((1+`r')^(11)) +(car_tax)/((1+`r')^(12))+(car_tax)/((1+`r')^(13)) +(car_tax)/((1+`r')^(14)) +(car_tax)/((1+`r')^(15))+(car_tax)/((1+`r')^(16))+(car_tax)/((1+`r')^(17))+(car_tax)/((1+`r')^(18))+(car_tax)/((1+`r')^(19))+(car_tax)/((1+`r')^(20))+(car_tax)/((1+`r')^(21))+(car_tax)/((1+`r')^(22))+(car_tax)/((1+`r')^(23))+(car_tax)/((1+`r')^(24))  if drivetype!=3 & env_cat>2

replace pv_tax_long=car_tax+(car_tax)/((1+`r')^1)+(car_tax)/((1+`r')^2)+(car_tax)/((1+`r')^3)+(car_tax/0.4)/((1+`r')^(4)) +(car_tax/0.4)/((1+`r')^(5)) +(car_tax/0.4)/((1+`r')^(6)) +(car_tax/0.4)/((1+`r')^(7)) +(car_tax/0.4)/((1+`r')^(8)) +(car_tax/0.4)/((1+`r')^(9)) +(car_tax/0.4)/((1+`r')^(10)) +(car_tax/0.4)/((1+`r')^(11)) +(car_tax/0.4)/((1+`r')^(12))+(car_tax/0.4)/((1+`r')^(13)) +(car_tax/0.4)/((1+`r')^(14))+(car_tax/0.4)/((1+`r')^(15))+(car_tax/0.4)/((1+`r')^(16))+(car_tax/0.4)/((1+`r')^(17))+(car_tax/0.4)/((1+`r')^(18))+(car_tax/0.4)/((1+`r')^(19))+(car_tax/0.4)/((1+`r')^(20))+(car_tax/0.4)/((1+`r')^(21))+(car_tax/0.4)/((1+`r')^(22))+(car_tax/0.4)/((1+`r')^(23))+(car_tax/0.4)/((1+`r')^(24)) if drivetype==3

replace pv_tax_long=car_tax+(car_tax)/((1+`r')^1)+(car_tax)/((1+`r')^2)+(car_tax)/((1+`r')^3)+(car_tax/0.6)/((1+`r')^(4)) +(car_tax/0.6)/((1+`r')^(5)) +(car_tax/0.6)/((1+`r')^(6)) +(car_tax/0.6)/((1+`r')^(7)) +(car_tax/0.6)/((1+`r')^(8)) +(car_tax/0.6)/((1+`r')^(9)) +(car_tax/0.6)/((1+`r')^(10)) +(car_tax/0.6)/((1+`r')^(11)) +(car_tax/0.6)/((1+`r')^(12))+(car_tax/0.6)/((1+`r')^(13)) +(car_tax/0.6)/((1+`r')^(14))+(car_tax/0.6)/((1+`r')^(15))+(car_tax/0.6)/((1+`r')^(16))+(car_tax/0.6)/((1+`r')^(17))+(car_tax/0.6)/((1+`r')^(18))+(car_tax/0.6)/((1+`r')^(19))+(car_tax/0.6)/((1+`r')^(20))+(car_tax/0.6)/((1+`r')^(21))+(car_tax/0.6)/((1+`r')^(22))+(car_tax/0.6)/((1+`r')^(23))+(car_tax/0.6)/((1+`r')^(24)) if drivetype!=3 & env_cat==1

replace pv_tax_long=car_tax+(car_tax)/((1+`r')^1)+(car_tax)/((1+`r')^2)+(car_tax)/((1+`r')^3)+(car_tax/0.8)/((1+`r')^(4)) +(car_tax/0.8)/((1+`r')^(5)) +(car_tax/0.8)/((1+`r')^(6)) +(car_tax/0.8)/((1+`r')^(7)) +(car_tax/0.8)/((1+`r')^(8)) +(car_tax/0.8)/((1+`r')^(9)) +(car_tax/0.8)/((1+`r')^(10)) +(car_tax/0.8)/((1+`r')^(11)) +(car_tax/0.8)/((1+`r')^(12))+(car_tax/0.8)/((1+`r')^(13)) +(car_tax/0.8)/((1+`r')^(14)) +(car_tax/0.8)/((1+`r')^(15))+(car_tax/0.8)/((1+`r')^(16))+(car_tax/0.8)/((1+`r')^(17))+(car_tax/0.8)/((1+`r')^(18))+(car_tax/0.8)/((1+`r')^(19))+(car_tax/0.8)/((1+`r')^(20))+(car_tax/0.8)/((1+`r')^(21))+(car_tax/0.8)/((1+`r')^(22))+(car_tax/0.8)/((1+`r')^(23))+(car_tax/0.8)/((1+`r')^(24)) if drivetype!=3 & env_cat==2
*Also generate the pv of driving costs: 
gen pv_drive_cost_long=drive_cost_ann + drive_cost_ann*((1-((1+`r')^(-(24))))/`r') 

gen var_cost_long=(pv_tax_long+pv_drive_cost_long)/1000

drop pv_tax_long pv_drive_cost_long

**Assuming a car-lifetime of 20 years and an interest rate of 6% for the sensitivity check mid_long:
local r=0.06

gen pv_tax_long=.

replace pv_tax_long=car_tax+(car_tax)/((1+`r')^1)+(car_tax)/((1+`r')^2)+(car_tax)/((1+`r')^3)+(car_tax)/((1+`r')^(4)) +(car_tax)/((1+`r')^(5)) +(car_tax)/((1+`r')^(6)) +(car_tax)/((1+`r')^(7)) +(car_tax)/((1+`r')^(8)) +(car_tax)/((1+`r')^(9)) +(car_tax)/((1+`r')^(10)) +(car_tax)/((1+`r')^(11)) +(car_tax)/((1+`r')^(12))+(car_tax)/((1+`r')^(13)) +(car_tax)/((1+`r')^(14)) +(car_tax)/((1+`r')^(15))+(car_tax)/((1+`r')^(16))+(car_tax)/((1+`r')^(17))+(car_tax)/((1+`r')^(18))+(car_tax)/((1+`r')^(19)) if drivetype!=3 & env_cat>2

replace pv_tax_long=car_tax+(car_tax)/((1+`r')^1)+(car_tax)/((1+`r')^2)+(car_tax)/((1+`r')^3)+(car_tax/0.4)/((1+`r')^(4)) +(car_tax/0.4)/((1+`r')^(5)) +(car_tax/0.4)/((1+`r')^(6)) +(car_tax/0.4)/((1+`r')^(7)) +(car_tax/0.4)/((1+`r')^(8)) +(car_tax/0.4)/((1+`r')^(9)) +(car_tax/0.4)/((1+`r')^(10)) +(car_tax/0.4)/((1+`r')^(11)) +(car_tax/0.4)/((1+`r')^(12))+(car_tax/0.4)/((1+`r')^(13)) +(car_tax/0.4)/((1+`r')^(14))+(car_tax/0.4)/((1+`r')^(15))+(car_tax/0.4)/((1+`r')^(16))+(car_tax/0.4)/((1+`r')^(17))+(car_tax/0.4)/((1+`r')^(18))+(car_tax/0.4)/((1+`r')^(19)) if drivetype==3

replace pv_tax_long=car_tax+(car_tax)/((1+`r')^1)+(car_tax)/((1+`r')^2)+(car_tax)/((1+`r')^3)+(car_tax/0.6)/((1+`r')^(4)) +(car_tax/0.6)/((1+`r')^(5)) +(car_tax/0.6)/((1+`r')^(6)) +(car_tax/0.6)/((1+`r')^(7)) +(car_tax/0.6)/((1+`r')^(8)) +(car_tax/0.6)/((1+`r')^(9)) +(car_tax/0.6)/((1+`r')^(10)) +(car_tax/0.6)/((1+`r')^(11)) +(car_tax/0.6)/((1+`r')^(12))+(car_tax/0.6)/((1+`r')^(13)) +(car_tax/0.6)/((1+`r')^(14))+(car_tax/0.6)/((1+`r')^(15))+(car_tax/0.6)/((1+`r')^(16))+(car_tax/0.6)/((1+`r')^(17))+(car_tax/0.6)/((1+`r')^(18))+(car_tax/0.6)/((1+`r')^(19))  if drivetype!=3 & env_cat==1

replace pv_tax_long=car_tax+(car_tax)/((1+`r')^1)+(car_tax)/((1+`r')^2)+(car_tax)/((1+`r')^3)+(car_tax/0.8)/((1+`r')^(4)) +(car_tax/0.8)/((1+`r')^(5)) +(car_tax/0.8)/((1+`r')^(6)) +(car_tax/0.8)/((1+`r')^(7)) +(car_tax/0.8)/((1+`r')^(8)) +(car_tax/0.8)/((1+`r')^(9)) +(car_tax/0.8)/((1+`r')^(10)) +(car_tax/0.8)/((1+`r')^(11)) +(car_tax/0.8)/((1+`r')^(12))+(car_tax/0.8)/((1+`r')^(13)) +(car_tax/0.8)/((1+`r')^(14)) +(car_tax/0.8)/((1+`r')^(15))+(car_tax/0.8)/((1+`r')^(16))+(car_tax/0.8)/((1+`r')^(17))+(car_tax/0.8)/((1+`r')^(18))+(car_tax/0.8)/((1+`r')^(19)) if drivetype!=3 & env_cat==2
*Also generate the pv of driving costs: 
gen pv_drive_cost_long=drive_cost_ann + drive_cost_ann*((1-((1+`r')^(-(19))))/`r') 

gen var_cost_midlong=(pv_tax_long+pv_drive_cost_long)/1000

drop pv_tax_long pv_drive_cost_long


**Assuming a car-lifetime of 6 years and an interest rate of 6% for the sensitivity check short:
local r=0.06

gen pv_tax_short=.

replace pv_tax_short=car_tax+(car_tax)/((1+`r')^1)+(car_tax)/((1+`r')^2)+(car_tax)/((1+`r')^3)+(car_tax)/((1+`r')^(4)) +(car_tax)/((1+`r')^(5)) if drivetype!=3 & env_cat>2

replace pv_tax_short=car_tax+(car_tax)/((1+`r')^1)+(car_tax)/((1+`r')^2)+(car_tax)/((1+`r')^3)+(car_tax/0.4)/((1+`r')^(4)) +(car_tax/0.4)/((1+`r')^(5))  if drivetype==3

replace pv_tax_short=car_tax+(car_tax)/((1+`r')^1)+(car_tax)/((1+`r')^2)+(car_tax)/((1+`r')^3)+(car_tax/0.6)/((1+`r')^(4)) +(car_tax/0.6)/((1+`r')^(5))   if drivetype!=3 & env_cat==1

replace pv_tax_short=car_tax+(car_tax)/((1+`r')^1)+(car_tax)/((1+`r')^2)+(car_tax)/((1+`r')^3)+(car_tax/0.8)/((1+`r')^(4)) +(car_tax/0.8)/((1+`r')^(5))  if drivetype!=3 & env_cat==2

*Also generate the pv of driving costs: 
gen pv_drive_cost_short=drive_cost_ann + drive_cost_ann*((1-((1+`r')^(-(5))))/`r') 

gen var_cost_short=(pv_tax_short+pv_drive_cost_short)/1000

drop pv_tax_short pv_drive_cost_short

**Assuming a car-lifetime of 10 years and an interest rate of 6% for the sensitivity check short:
local r=0.06

gen pv_tax_short=.

replace pv_tax_short=car_tax+(car_tax)/((1+`r')^1)+(car_tax)/((1+`r')^2)+(car_tax)/((1+`r')^3)+(car_tax)/((1+`r')^(4)) +(car_tax)/((1+`r')^(5))+(car_tax)/((1+`r')^(6))+(car_tax)/((1+`r')^(7))+(car_tax)/((1+`r')^(8)) +(car_tax)/((1+`r')^(9)) if drivetype!=3 & env_cat>2

replace pv_tax_short=car_tax+(car_tax)/((1+`r')^1)+(car_tax)/((1+`r')^2)+(car_tax)/((1+`r')^3)+(car_tax/0.4)/((1+`r')^(4)) +(car_tax/0.4)/((1+`r')^(5)) +(car_tax/0.4)/((1+`r')^(6))+(car_tax/0.4)/((1+`r')^(7))+(car_tax/0.4)/((1+`r')^(8))+(car_tax/0.4)/((1+`r')^(9)) if drivetype==3

replace pv_tax_short=car_tax+(car_tax)/((1+`r')^1)+(car_tax)/((1+`r')^2)+(car_tax)/((1+`r')^3)+(car_tax/0.6)/((1+`r')^(4)) +(car_tax/0.6)/((1+`r')^(5)) +(car_tax/0.6)/((1+`r')^(6))+(car_tax/0.6)/((1+`r')^(7))+(car_tax/0.6)/((1+`r')^(8))+(car_tax/0.6)/((1+`r')^(9))   if drivetype!=3 & env_cat==1

replace pv_tax_short=car_tax+(car_tax)/((1+`r')^1)+(car_tax)/((1+`r')^2)+(car_tax)/((1+`r')^3)+(car_tax/0.8)/((1+`r')^(4)) +(car_tax/0.8)/((1+`r')^(5)) +(car_tax/0.8)/((1+`r')^(6)) +(car_tax/0.8)/((1+`r')^(7)) +(car_tax/0.8)/((1+`r')^(8)) +(car_tax/0.8)/((1+`r')^(9))  if drivetype!=3 & env_cat==2

*Also generate the pv of driving costs: 
gen pv_drive_cost_short=drive_cost_ann + drive_cost_ann*((1-((1+`r')^(-(5))))/`r')

gen var_cost_midshort=(pv_tax_short+pv_drive_cost_short)/1000

drop pv_tax_short pv_drive_cost_short




reg price height weight var_cost typkw ev hyrid diesel env_friendly car_size brand_FE2 brand_FE3 brand_FE4 brand_FE5 brand_FE6 brand_FE7 brand_FE8 brand_FE9 brand_FE10 brand_FE11 brand_FE12 brand_FE13 brand_FE14 brand_FE15 brand_FE16 car_cat_FE*  dist_car_size dist_car_height dist_TYPKW dist_weight  hourly_wage
predict resid, residual

estimates store cont1
estadd local brand_fe "Yes"
estadd local type_fe "Yes"

test dist_car_size dist_car_height dist_TYPKW dist_weight  hourly_wage  
estadd local f_stat=round(r(F),.001)

reg price height weight var_cost typkw ev hyrid diesel car_size env_friendly brand_FE2 brand_FE3 brand_FE4 brand_FE5 brand_FE6 brand_FE7 brand_FE8 brand_FE9 brand_FE10 brand_FE11 brand_FE12 brand_FE13 brand_FE14 brand_FE15 brand_FE16 car_cat_FE*   dist_ind_car_size dist_ind_car_height dist_ind_TYPKW dist_ind_weight hourly_wage

predict resid_inddist, residual

estimates store cont2
estadd local brand_fe "Yes"
estadd local type_fe "Yes"

test dist_ind_car_size dist_ind_car_height dist_ind_TYPKW dist_ind_weight hourly_wage

estadd local f_stat=round(r(F),.001)


reg price height weight var_cost typkw ev hyrid diesel car_size env_friendly brand_FE2 brand_FE3 brand_FE4 brand_FE5 brand_FE6 brand_FE7 brand_FE8 brand_FE9 brand_FE10 brand_FE11 brand_FE12 brand_FE13 brand_FE14 brand_FE15 brand_FE16 car_cat_FE*   dist_car_size dist_car_height dist_TYPKW dist_weight  

predict resid_dist2, residual
estimates store cont3
estadd local brand_fe "Yes"
estadd local type_fe "Yes"

test dist_car_size dist_car_height dist_TYPKW dist_weight 

estadd local f_stat=round(r(F),.001)


reg price height weight var_cost_low typkw ev hyrid diesel car_size env_friendly brand_FE2 brand_FE3 brand_FE4 brand_FE5 brand_FE6 brand_FE7 brand_FE8 brand_FE9 brand_FE10 brand_FE11 brand_FE12 brand_FE13 brand_FE14 brand_FE15 brand_FE16 car_cat_FE*   dist_car_size dist_car_height dist_TYPKW dist_weight  hourly_wage 

predict resid_low, residual
estimates store cont4
estadd local brand_fe "Yes"
estadd local type_fe "Yes"

test dist_car_size dist_car_height dist_TYPKW dist_weight  hourly_wage   
estadd local f_stat=round(r(F),.001)

reg price height weight var_cost_high typkw ev hyrid diesel car_size env_friendly brand_FE2 brand_FE3 brand_FE4 brand_FE5 brand_FE6 brand_FE7 brand_FE8 brand_FE9 brand_FE10 brand_FE11 brand_FE12 brand_FE13 brand_FE14 brand_FE15 brand_FE16 car_cat_FE*   dist_car_size dist_car_height dist_TYPKW dist_weight  hourly_wage 

predict resid_high, residual
estimates store cont5
estadd local brand_fe "Yes"
estadd local type_fe "Yes"

test dist_car_size dist_car_height dist_TYPKW dist_weight  hourly_wage  
estadd local f_stat=round(r(F),.001)


reg price height weight var_cost_short typkw ev hyrid diesel car_size env_friendly brand_FE2 brand_FE3 brand_FE4 brand_FE5 brand_FE6 brand_FE7 brand_FE8 brand_FE9 brand_FE10 brand_FE11 brand_FE12 brand_FE13 brand_FE14 brand_FE15 brand_FE16 car_cat_FE*   dist_car_size dist_car_height dist_TYPKW dist_weight  hourly_wage 

predict resid_short, residual
estimates store cont6
estadd local brand_fe "Yes"
estadd local type_fe "Yes"

test dist_car_size dist_car_height dist_TYPKW dist_weight  hourly_wage  
estadd local f_stat=round(r(F),.001)

reg price height weight var_cost_midshort typkw ev hyrid diesel car_size env_friendly brand_FE2 brand_FE3 brand_FE4 brand_FE5 brand_FE6 brand_FE7 brand_FE8 brand_FE9 brand_FE10 brand_FE11 brand_FE12 brand_FE13 brand_FE14 brand_FE15 brand_FE16 car_cat_FE*   dist_car_size dist_car_height dist_TYPKW dist_weight  hourly_wage 

predict resid_midshort, residual
estimates store cont6_1
estadd local brand_fe "Yes"
estadd local type_fe "Yes"

test dist_car_size dist_car_height dist_TYPKW dist_weight  hourly_wage  
estadd local f_stat=round(r(F),.001)

reg price height weight var_cost_long typkw ev hyrid diesel car_size env_friendly brand_FE2 brand_FE3 brand_FE4 brand_FE5 brand_FE6 brand_FE7 brand_FE8 brand_FE9 brand_FE10 brand_FE11 brand_FE12 brand_FE13 brand_FE14 brand_FE15 brand_FE16 car_cat_FE*   dist_car_size dist_car_height dist_TYPKW dist_weight  hourly_wage 

predict resid_long, residual
estimates store cont7
estadd local brand_fe "Yes"
estadd local type_fe "Yes"

test dist_car_size dist_car_height dist_TYPKW dist_weight  hourly_wage  
estadd local f_stat=round(r(F),.001)

reg price height weight var_cost_midlong typkw ev hyrid diesel car_size env_friendly brand_FE2 brand_FE3 brand_FE4 brand_FE5 brand_FE6 brand_FE7 brand_FE8 brand_FE9 brand_FE10 brand_FE11 brand_FE12 brand_FE13 brand_FE14 brand_FE15 brand_FE16 car_cat_FE*   dist_car_size dist_car_height dist_TYPKW dist_weight  hourly_wage 

predict resid_midlong, residual
estimates store cont7_1
estadd local brand_fe "Yes"
estadd local type_fe "Yes"

test dist_car_size dist_car_height dist_TYPKW dist_weight  hourly_wage  
estadd local f_stat=round(r(F),.001)


reg price height weight var_cost_variable typkw ev hyrid diesel car_size env_friendly brand_FE2 brand_FE3 brand_FE4 brand_FE5 brand_FE6 brand_FE7 brand_FE8 brand_FE9 brand_FE10 brand_FE11 brand_FE12 brand_FE13 brand_FE14 brand_FE15 brand_FE16 car_cat_FE*  dist_car_size dist_car_height dist_TYPKW dist_weight  hourly_wage 

predict resid_variable, residual
estimates store cont8
estadd local brand_fe "Yes"
estadd local type_fe "Yes"

test dist_car_size dist_car_height dist_TYPKW dist_weight  hourly_wage  
estadd local f_stat=round(r(F),.001)



#delimit;
esttab cont1 cont2 cont3 cont4 cont5 cont6 cont6_1 cont7 cont7_1 cont8 using "$root/results/additional_results/control_function03_complete.tex", booktabs replace  rename(var_cost_low var_cost var_cost_variable var_cost var_cost_midlong var_cost var_cost_midshort var_cost var_cost_long var_cost var_cost_short var_cost var_cost_high var_cost) keep(dist_car_size dist_car_height dist_TYPKW dist_weight  hourly_wage dist_ind_car_size dist_ind_car_height dist_ind_TYPKW dist_ind_weight height weight  typkw ev hyrid diesel car_size env_friendly var_cost) 
mgroups( "Base" "IV different" "Var cost robustness", pattern(1 1 0 1 0 0 0 0 0 0 0) span lhs(Specification: ) prefix(\multicolumn{@span}{c}{) suffix(})erepeat(\cmidrule(lr){@span})) cells(b(star fmt(3)) se(fmt(2) par)) collabels(none)  starlevels(+ 0.1 * 0.05 ** 0.01 *** 0.001) 
stats(N r2 aic brand_fe type_fe f_stat, fmt(%9.0g)  label("Observations" "$R^2$" "AIC" "Brand dummy" "Car type dummy" "F-stat (relevance)"))  order(dist_car_size dist_car_height dist_TYPKW dist_weight  hourly_wage dist_ind_car_size dist_ind_car_height dist_ind_TYPKW dist_ind_weight height weight  typkw ev hyrid diesel car_size env_friendly var_cost)
scalars() coeflabels(typkw "Engine power (KW)" height "Car height" weight "Car weight" car_size "Car size" var_cost "Variable costs" hourly_wage "Wage manufacturing"  diesel "Diesel engine" hyrid "Hybrid engine" ev "Electric engine" env_friendly "Below category C " dist_car_size "Distance size" dist_ind_car_size "Distance (Nb.) car size" dist_TYPKW "Distance Engine power" dist_ind_TYPKW "Distance (Nb.) engine power" dist_weight "Distance weight" dist_ind_weight "Distance (Nb.) weight" dist_car_height "Distance height" dist_ind_car_height "Distance (Nb.) height"  )  
alignment(D{.}{.}{-1}) 
;

#delimit cr



#delimit;
esttab cont1 cont2 cont3 cont4 cont7 cont8 using "$root/results/tables_graphs/appendixcontrol_function03.tex", booktabs replace  rename(var_cost_low var_cost var_cost_variable var_cost var_cost_long var_cost) keep(dist_car_size dist_car_height dist_TYPKW dist_weight  hourly_wage dist_ind_car_size dist_ind_car_height dist_ind_TYPKW dist_ind_weight height weight  typkw ev hyrid diesel car_size env_friendly var_cost) 
mgroups( "Base" "IV different" "Var cost robustness", pattern(1 1 0 1 0 0) span lhs(Specification: ) prefix(\multicolumn{@span}{c}{) suffix(})erepeat(\cmidrule(lr){@span})) cells(b(star fmt(3)) se(fmt(2) par)) collabels(none)  starlevels(+ 0.1 * 0.05 ** 0.01 *** 0.001) stats(N r2 aic brand_fe type_fe f_stat, fmt(%9.0g)  label("Observations" "$R^2$" "AIC" "Brand dummy" "Car type dummy" "F-stat (relevance)"))  order(dist_car_size dist_car_height dist_TYPKW dist_weight  hourly_wage dist_ind_car_size dist_ind_car_height dist_ind_TYPKW dist_ind_weight height weight  typkw ev hyrid diesel car_size env_friendly var_cost)
scalars() coeflabels(typkw "Engine power (KW)" height "Car height" weight "Car weight" car_size "Car size"  diesel "Diesel engine" hyrid "Hybrid engine" ev "Electric engine" env_friendly "Below Category C " var_cost "Variable costs" hourly_wage "Wage manufacturing" dist_car_size "Distance size" dist_ind_car_size "Distance (Nb.) car size" dist_TYPKW "Distance Engine power" dist_ind_TYPKW "Distance (Nb.) engine power" dist_weight "Distance weight" dist_ind_weight "Distance (Nb.) weight" dist_car_height "Distance height" dist_ind_car_height "Distance (Nb.) height" )  
alignment(D{.}{.}{-1}) 
;

#delimit cr



**Drop the no longer needed variables (instruments) to have a smaller dataset in memory: 


drop MAKName hourly_wage dist_car_size dist_ind_car_size dist_car_height dist_ind_car_height dist_TYPKW dist_ind_TYPKW dist_weight dist_ind_weight dist_price dist_ind_price  pv_drive_cost pv_drive_cost_var _est_cont1 _est_cont2 _est_cont3 _est_cont4 _est_cont5 _est_cont6 _est_cont7 _est_cont6_1 _est_cont7_1 _est_cont8


*Generate an id variable to runs from 1 to N: 
egen id=group(PID)
*Generate a car_option variable that runs from 1 to J: 
egen car_option=group(car_combo)

*Drop if there are any missing values, otherwise matlab script does not work - the missing values are only in the car variables, as not all options were available in every year - hence we drop based on price:
drop if missing(price)

*Order the variables as used in the estimation: 
order id car_option choice price height weight var_cost typkw ev hyrid diesel car_size size_hh2 size_hh3 size_hh4 size_hh5 age_kw2 age_kw3 ev_urban2 ev_urban3 ev_dist_ev ev_range5 ev_year2 ev_year3 ev_pv ev_ho ev_wealth2 ev_wealth3 ev_wealth4 price_wealth2 price_wealth3 price_wealth4 brand_FE2 brand_FE3 brand_FE4 brand_FE5 brand_FE6 brand_FE7 brand_FE8 brand_FE9 brand_FE10 brand_FE11 brand_FE12 brand_FE13 brand_FE14 brand_FE15 brand_FE16 car_cat_FE2 car_cat_FE3 car_cat_FE4 car_cat_FE5 car_cat_FE6 car_cat_FE7 car_cat_FE8 env_friendly 

*Order the additionally variables resid and var_costs which are used for robustness checks:
order resid*, after(env_friendly)
order var_cost_*,after(resid_variable)


*Save the file as csv to estimate the models in Matlab:

preserve 

*drop the all non-required vars such that memory is as small as possible: 
drop year_reg PID car_combo drivetype et_verbrauch et_co2 env_cat prhnewprice average_use_var drive_cost drive_cost_ann_variable car_tax TOT_EINK TOT_EINK_orig TOTVERM wealth wealth_group el_verbrauch drive_cost_ann wealth_quart2 pv_tax wealth_quart2 wealth_quart3 wealth_quart4 average_use


gen choice_sit=id 

order choice_sit, after(id)
order choice, after(choice_sit)

sort id car_option


export delimited using "$root/Code/Estimation/matlab/EV_project09.csv", replace novarnames

restore

**Create the welfare dataset to use once we have the estimated parameters: 

fillin id car_option

preserve

tempfile socios

gcollapse (firstnm) PID year_reg TOT_EINK TOT_EINK_orig TOTVERM wealth_group wealth, by(id)

save `socios' 

restore 


merge m:1 id using `socios', update


drop _merge resid_* var_cost_*  _fillin


*Create the wealth_group: 
tab wealth_group, g(wealth_quart)
drop wealth_quart1



**Save the dataset as EV_Welfare04.dta, which is the dataset that is used to produce all the post-regression estimates:

save "$root/Data/Produced/EV_welfare04.dta", replace



